Anharmonic phonon frequency shift in MgB2 
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We compute the anharmonic shift of the phonon frequencies in MgB2, using density functional 
theory. We explicitly take into account the scattering between different phonon modes at different 
q-points in the Brillouin zone. The shift of the E2 9 mode at the F point is +5% of the harmonic 
frequency. This result comes from the cancellation between the contributions of the four- and 
three-phonon scattering, respectively +10% and —5%. A similar shift is predicted at the A point, 
in agreement with inelastic X-ray scattering phonon-dispersion measurements. A smaller shift is 
observed at the M point. 

PACS numbers: 63.20.Dj, 63.20.Kr, 71. 15. Mb 



The discovery of high critical temperature (T c = 
39K) in MgB 2 has challenged our understanding 
of electron-phonon coupling mediated superconductivity. 
Several mechanisms might cooperate in determining the 
actual value of T c in MgBo, including the existence of a 
double gap structure 0, 0, , and anharmonic effects 
[M IE M EJ IHJ • Much attention has been devoted to 
the study of the E 2g phonon mode, which consists of an 
antiphase vibration of the two boron atoms, parallel to 
the hexagon plane. This mode has the largest coupling 
with electrons according to inelastic X-ray scattering 
measurements 11], in agreement with a commonly ac- 
cepted interpretation of the MgB 2 electronic band struc- 
ture Il2t Il3|. The E 2g mode has also been associated 
with a supposedly large anharmonicity, but a quantita- 
tive determinations of the anharmonic effects has led so 
far to controversial results @, H S H H E3 ■ 

In particular, actual theoretical calculations indicate 
the anharmonic frequency shift of the Ei g mode at T 
to be large, with varying estimates that range from 
+15% HQ, up to +20/25% HE1 of the theoretical har- 
monic frequency (~ 65 meV). A comparably large shift 
is expected to affect this mode all al ong the T-A direc- 
tion Raman measurements [I3 . ll5Lll6lll^ | seem to 
confirm the prediction of a lar ge a nharmonicity at T. In 
fact, the best resolved spectra [lj|, present a peak at 77 
meV (i.e. 18% higher than the harmonic frequency) that 
could be attributed to the E 2g mode at T. However, the 
structure observed below 40 meV 0], and the clearly 
asymmetric profile (reminiscent of a Fano resonance) of 
the 77 meV peak indicate that the Raman experiment 
probes not just phonon vibrations, but electronic excita- 
tions as well. Thus, the position of the 77 meV peak does 
not necessarily correspond to the E 2g phonon frequency. 
Moreover, Raman spectra display a drastic dependence 
on the temperature. The width of the 77 meV peak is 
~ 20 meV at T= 40K and reaches almost 40 meV at room 
temperature 0, 0] . This behavior has been tentatively 
attributed to anharmonicity [lfil Il7j. but, according to 
the calculations of Ref . [Jjj , the E 2g -F anharmonic broad- 
ening at room temperature is negligible, just 1.2 meV. 



The determination of the phonon frequency from Raman 
data requires further analysis, as suggested in Ref. [ljj, 
together with a clear understanding of the temperature 
dependence. 

On the other hand, the prediction of a strong anhar- 
monic frequency shift is in stark contrast with the inelas- 
tic X-ray measurements of the MgB 2 phonon dispersion 
of Ref. [llj , which E 2g phonon frequencies measurements 
are in agreement (within 5%) with the calculated har- 
monic phonon frequencies 0, 0, Q] at several points 
along the T-A direction. Contrary to Raman scatter- 
ing |l6llr7j|. in X-ray inelastic scattering the cross-section 
of electronic excitations is much smaller than that of the 
phonon excitations. Thus, X-ray measurements are a di- 
rect probe of the phonon frequencies. Small anharmonic 
frequency shift is also suggested by a recent point-contact 
spectroscopy measurement [Tij ]. 

The origin of the discrepancy between thepredicted 
large anharmonic shift of the E 2g mode |5j, |g, [jj, |a 1 9,11(1 
and the inelastic X-ray scattering measurements |llj . 
might be traced back to the approximations involved in 
the calculations. We recall here that the anharmonic fre- 
quency shift of a phonon mode (qj) (identified by the 
reciprocal-space vector q and by the band index j) is 
due to the anharmonic interaction with the ensemble of 
all phonons, having different momentum q and different 
j. Up to now, theoretical studies of MgB 2 have been 
based on the study of a frozen-phonon energy profile. 
This approach completely neglects: (i) the dependence of 
the anharmonic scattering on the exchanged phonon mo- 
mentum q, and (ii) the coupling between phonons having 
different j. The calculation of the anharmonic scatter- 
ing between phonons having different q and different j 
is possible, using first-principles [13] techniques. E.g., 
it has allowed allowed the determination of anharmonic 
line shift and broadening of the Raman modes in covalent 
semiconductors (lif . However, this approach is compu- 
tationally demanding and has not been yet applied to 
MgB 2 . 

In this Letter we compute the anharmonic frequency 
shift of the phonon E 2g mode of MgB 2 from first- 
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FIG. 1: Diagrammatic representation of the leading terms in 
the perturbative expansion of a non-harmonic phonon hamil- 
tonian (see the text). 



principles |20|. We explicitly take into account the scat- 
tering between different phonon modes at different q- 
points in the reciprocal space. We use the density- 
functional perturbation theory of Ref. [2l| ■ The anhar- 
monic terms of the phonon Hamiltonian are obtained us- 
ing a method, based on the "2n + 1" theorem, which 
allows the efficient calculation of third order derivatives 
of total energy, at any point in the reciprocal space, in 
metallic systems [iij . 

The phonon excitations of a crystal are determined by 
the interatomic energy, and correspond to the poles of the 
interacting phonon Green function [2^, 01- n qj (w) rep- 
resents the anharmonic contribution to the phonon self- 
energy of the phonon mode (qj) having frequency Lo q j. 
In the case \H q j\ <C u> q j, the frequency shift A q j and line 
broadening T q j, can be obtained as A q j = lZ e [H q j(uj q j)] 



and T q j = —2 m [Tl q j(uj q j)], lZ e andX m being the real and 
imaginary part of a complex number. The three lowest 
order terms in the perturbative expansion, for a fixed cell 
geometry, are: 



n 4j ( w ) = ifi'( u )+C( w ) + n^( w ) I 



r(-B) 



(1) 



which correspond, respectively, to the tadpole (T), loop 
(L), and bubble (B) diagrams of Fig. whose calculation 
is described in e.g. in Refs. |23ll24f. 

Let us consider the interatomic potential energy 
£ tot ({u sa (Ri)}), where u sa (Ri) is the displacement from 
the equilibrium position of the s-th atom in the crystal 
cell identified by the lattice vector B, along the a Carte- 
sian coordinate. We define 



n ccell 



d n £ 



where £ cel1 is the energy per unit cell and the 
adimcnsional quantity u q j is defined by u q j — 

-T 



2M. 



■v sa (qj)u sa {Ri)e lq R ', v sa (qj) being 



the orthogonal phonon eigenmodes normalized on the 
unit cell, M s the atomic mass, and N the number of 
q-points describing the system (or unit cells). Using this 
definition y( n ) is a intensive quantity, independent of N, 
at any order n. The evaluation of the three diagrams for 
the j zone center mode (q = 0) yields: 



q, 31,32 

n (s)/, a 1 \' i T/ (3)/ n - ■ „ ■ m2 ( gHji + ^X 1 + ngh + n qh ) 2(uj qjl - uj qj2 )(n qh - n qjl ) 



(2) 



where is a sum on the Brillouin zone (BZ), n q j is the 
Bose-statistics occupation of the phonon mode (qj), and 
S is a small positive number. Interpreting the three dia- 
grams of Fig^fii terms of scattering between phonons, 
the self-energy IL q j(u) can thus be expressed as a sum, 
over the BZ, of the phonon scattering coefficients V^ n \ 

Eqs. 1 1121 are valid at fixed cell geometry. In general, 
the dependence of YL q j(uj) on the temperature depends 
on further terms related to the thermal expansion of the 
lattice. However, if one can determine the lattice ther- 
mal expansion independently (e.g. experimentally), the 
shift at a given temperature can be obtained computing 



the harmonic frequency at the lattice parameters corre- 
sponding to that temperature. While only the term 
in Eq. ^has an imaginary component (i.e. is providing 
a contribution to the line broadening, which we already 
calculated and discussed in Ref. 11]), all three diagrams 
have real part and contribute to the shift. Due sym- 
metry, the Il( T ) term is zero j2^| in MgB2, and the fre- 
quency shift can, thus, be decomposed in three- and four- 

(3) ( B} 

phonon scattering contributions Aq ■ — lZ e \Tv Qj ■ (uoj)] 
and A^=H e [n^(co oj )}. 

Several authors used the frozen-phonon approach to 
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A (4) 


A (4) 


A (3) 


tot 


9.46 

(+15%) 


T = K, 
6.00 
(+9.2%) 


u = 65.14 meV 
-2.88 
(-4.4%) 


3.12 

(+4.8%) 


8.91 

(+14%) 


7.18 
(+11%) 


-2.99 
(-4.6%) 


4.19 
(+6.4%) 


10.68 

(+16%) 


T = 300 K 
6.52 
(+10%) 


lu = 64.51 meV 
-3.03 
(-4.7%) 


3.50 
(+5.4%) 


10.01 

(+15%) 


7.77 

(+12%) 


-3.12 

(-4.8%) 


4.65 
(+7.2%) 



TABLE I: Phonon frequency shift (meV) of the E 29 mode 
at the BZ center (T point). A (3) (A (4) ) is the shift due to 
the three(four)-phonon scattering contribution. A tot (A ° = 
A (3) + A (4) ) is the total shift. A$ D is defined in the text. 
uj is the unperturbed phonon frequency, and the values in 
parentheses are the shift as a percentage with respect to the 
u) at the corresponding temperature (T). Bold face values are 
obtained using the 16 x 16 x 16 electronic integration grid, 
the others from a less converged 14 x 14 x 8 grid [3l| . 

determine MgB 2 frequency shift H S H @- This 
approach is based on the study of the energy pro- 
file £g=o,j(A) obtained displacing MgB 2 atomic posi- 
tions by a value proportional to A, along the eigenvec- 
tor of mode j at q = 0. This approach is valid when: 
(i) the coupling of the mode j with different modes 
is negligible; (ii) the anharmonic interaction does not 
change over the BZ, i.e. when (Oj, q 2 j, q n j) ~ 
V( n \0j,0j,...,0j) Vq 2 ,...,q n eBZ. Only in this non- 
dispersive regime the behavior of £q=o,i(A) gives quanti- 
tative informations about the shift. 

Our electronic structure calculations is performed us- 
ing a plane-wave basis set up to a 35 Ry cutoff, and us- 
ing the pseudopotential approach j2(|. The terms V^ 3 ) 
of Eq . are computed analytically with the approach of 
Ref. 22] and the are obtained by finite differentia- 
tion of . Calculations are repeated both at T=0 K 
and T=300 K temperature crystal structures First 
order Hermite-Gaussian smearing of 0.025 Ry is used. 
The dynamical matrices and phonon frequencies are cal- 
culated using the PBE 29] approximation for the ex- 
change and correlation functional and a well converged 
16 x 16 x 16 Monkhorst-Pack grid for the electronic BZ 
integration. Due to the actual implementation in our 
code, third and fourth order interatomic force constants 
(V^ and 4) ) are obtained using the LDA 30] func- 
tional. By performing frozen-phonon calculations at T 
using both PBE and LDA we estimate that the total 
shift obtained by LDA should be slightly larger (< 0.5 
meV) than with PBE. 

In Tab. [I] we show the computed values of A^ 4 ) and 
A (3) (previously defined) for the E2 9 mode at T. To 





w 


A (4) 


A (3) 


j/\tot 


r,E 2g 

A (E 2fl ) 
M (E 2g ) 
M (E 29 ) 


65.14 
57.50 
87.18 
93.10 


T = K 

7.18 

7.77 

4.12 

1.10 


-2.99 
-3.10 
-1.98 
-1.29 


4.19 (+6.4%) 
4.67 (+8.1%) 
2.13 (+2.4%) 
-0.19 (-0.2%) 


r,Ei„ 
r,A 2u 
r,Bi s 


40.03 
46.21 
84.49 


0.48 
1.45 
-1.02 


-0.24 
-0.36 
-0.16 


0.25 (+0.6%) 
1.09 (+2.4%) 
-1.18 (-1.4%) 


r,E 29 

A (E 29 ) 
M (E 29 ) 
M (E 29 ) 


64.51 
56.88 
86.65 
92.80 


T = 300 K 
7.77 
8.16 
4.65 
1.20 


-3.12 
-3.50 
-2.13 
-1.46 


4.65 (+7.2%) 

4.66 (+8.2%) 
2.52 (+2.9%) 
-0.26 (-0.3%) 


r,Ei u 
r,A 2 „ 
r,Bi 9 


39.74 
45.81 
84.28 


0.76 
0.20 
-1.03 


-0.30 
-0.64 
-0.21 


0.46 (+1.2%) 
-0.45 (-1.0%) 
-1.24 (-1.5%) 



TABLE II: Phonon frequency shift (meV) of various modes 
at the three high symmetry points F, A, and M of the BZ 
(see also the caption of Tab. [IJ. The shifts of this table are 
obtained using the 14 x 14 x 8 electronic integration grid |3l|. 



make a comparison with the frozen-phonon approach, 
in Tab. [I] we also show A^^>, namely the shift due to 
the four-phonon scattering obtained neglecting the in- 
teraction with the other modes, and assuming a non- 
dispersive anharmonic potential: i. e. imposing in Eq. |3 
^(Oj.Oj.-gji.gji) = V^(0j,0j,0j,0j)5 ju j. As ex- 
pected 0, IE S S HD ! ^nd gi yes a large positive shift 
(+16% of the unperturbed frequency at T=300K), which, 
however, is far from being a quantitative estimate of the 
shift. In fact, the actual determination of scattering pro- 
cesses all over the BZ, and the inclusion of the other 
modes (ji ^ j), reduces the value of the four-phonon 
shift (AW in Tab. UJ) to just +10%. Furthermore, the 
inclusion of the A^ 3 ) term (which has a negative sign) 
reduces the total shift to just +5.4%. The main reason 
of the negative sign of A^ 3 ) is that the joint density of 
phonon states p(ui) — J2 q j 1 j 2 ^(^^^qji ~ u qj-i) i s almost 
entirely (more than 90%) distributed above the frequency 
of the r-E2 g mode. This implies that, in the summation 
done to obtain TZ e (U^) from Eq.EJat T= 0, the largest 
part of the terms have a negative value. 

To compute the shift at the A and M BZ points (q ^ 0) 
we use, respectively, a 1 x 1 x 2, and a 2 x 2 x 1 super-cell. 
Due to the heavier computational effort required, 
and are calculated with less converged electronic 

integrati on g rids equivalent to a 14 x 14 x 8 grid for the 1 x 
lxl cell |3l|. From Tab.^one can infer that calculations 
with this grid provide a qualitative description of the 
system. In Tab. ITT1 we show the shift, at A and M, of 
the modes belonging to the E 2g band, and, at T, the 
shift of modes different than E2 9 . It is apparent that 
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the total E 2g shifts at T and A are greater than both 
the E 2g shift at M, and the shift of the other modes at 
r. This is not surprising: it is related to the strong 
electron-phonon coupling of the E 2g mode along the T-A 
direction , and confirms one accepted interpretation of 
the MgB 2 electronic band structure 0, 0. Moreover, 
the E 2g shift at T and A are very similar, being 4.65 and 
4.66 meV (Tab. |Hj|, respectively, at room temperature. 
This confirms the predicted presence of a uniform shift 
along the T- A direction U • 

Keeping in mind that our best estimate for the T shift 
is just +3.5 meV (i.e. +5.4% of the harmonic value, 
from Tab. HJ we expect the anharmonic shift to be of 
the same order all along T-A. Such a value is compati- 
ble with the room temperature inelastic scattering X-ray 
measurements of Ref. 01 ■ The measured E 2g disper- 
sion along the T-A direction, from Ref. is in agree- 
ment with harmonic model calculations. Indeed, com- 
paring the experimental E 2g frequency at A (59.0 meV) 
with our calculation, even including anharmonic effects 
(56.88 + 4.66 = 61.5 meV, using the 14 x 14 x 8 grid) the 
agreement remains well within the experimental error- 
bar and the typical DFT error. 

Finally, we notice that our calculations predict an al- 
most negligible dependence of the E 2g -r frequency from 
the temperature (being 65.14 + 3.12 = 68.26 meV at 
T= OK, and 64.51 + 3.50 = 68.01 meV at T= 300K). 
In fact, the largest part of this shift is due to the cou- 
pling with phonons having a Debye temperature larger 
than 300K (to > 26 meV). 

In conclusion, we have obtained the anharmonic 
phonon frequency shift in MgB 2 using density functional 
theory. Two ingredients which turned out to be essen- 
tial for a quantitative description of the anharmonic E 2g - 
mode frequency-shift are: (i) the explicit calculation of 
the scattering processes all over the BZ; (ii) the inclu- 
sion in the perturbative expansion of both the three- and 
four-phonon scattering contributions, which have oppo- 
site sign. The resulting phonon frequency shift at T is 
just +5.4% of the harmonic frequency, and is expected 
to be of the same order along the T-A direction, in agree- 
ment with inelastic X-ray scattering [ill ]. 
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